NASA Technical Memorandum 84402 ^84“24276 


An Explicit Predictor - Corrector 
Solver with Applications to Burgers’ 
Equation 

Suhrit K. Dey and Charlie Dey 


September 1983 


NASA 

National Aeronautics and 
Space Administration 




NASA Technical Memorandum 84402 


An Explicit Predictor- Corrector 
Solver with Applications to Burgers’ 
Equation 

Suhrit K. Dey, Ames Research Center, Moffett Field, California 
Charlie Dey, Country Lane School, San Jose, California 


NASA 

National Aeronautics and 
Space Administration 

Ames Research Center 

Moffett Field. California 94035 




AN EXPLICIT PREDICTOR-CORRECTOR SOLVER WITH APPLICATION 


TO BURGERS ' EQUATION 
Suhrit K. Dey and Charlie Dey* 
Ames Research Center 


SUMMARY 


Forward Euler's explicit, finite-difference formula of extrapolation is used as 
a predictor and a convex formula as a corrector to integrate differential equations 
numerically. An application has been made to Burgers' equation. 


INTRODUCTION 


The single-step, explicit, forward Euler scheme is possibly the simplest algo- 
rithm with which to solve initial-value problems. However, it has vpry restricted 
stability properties which require small step sizes so that this formula of extrapo- 
lation may be useful. Implicit schemes, in which large step sizes may be used, have 
better stability properties. But iterative methods are generally used for computa- 
tional solution. Most predictor-corrector algorithms require multistep operations 
in which correctors use iterative methods for convergence. 

In this work, a filtering formula has been introduced by the second author so 
that extrapolated values obtained by an explicit formula may be appropriately cor- 
rected using one iteration which requires no computation of derivatives or inversion 
of a matrix. Mathematically, the corrector introduces a convex mapping, and as such 
it may be called a convex corrector. The explicit formula of extrapolation used in 
this work is the forward Euler difference scheme. Linearized stability analysis by 
Lomax (private communication) showed that the method is stable for step sizes that 
could be much larger than those required for stability of the forward Euler scheme. 
This will be discussed further on. With regard to the solution of nonlinear models, 
the limitations of this scheme are still being investigated. 

Let us first develop the algorithm, discuss its stability properties, and look 
into some of its simple applications. 


THE ALGORITHM 


Let us consider an initial-value model: 


du 

dt 


f (u,t) 


u(t ) 
o 


( 1 ) 


*Country Lane School, San Jose, Calif. 
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The forward Euler scheme is: 


( 2 ) 


U 


n+1 


= U + At f(U ,t ) 
n n n 


where At = time-step 

U = U(t ) = the net function corresponding to u(t ) 
n n n 

We may now set up a predictor-corrector formula as follows: 

U = U + At f (U ,t ) (predictor) (3) 

n n n 

U , = (1 - y)U + y[U + At f(U,t )] (corrector) (4) 

n+ 1 n n+ 1 

Gamma is herein called a filtering parameter. Here we have restricted y to 

0 < y < 1 (5) 

Thus, y is essentially a convex parameter. If y = 0, U n+1 = U. The corrector is 
ineffective. 

If the model (eq. (1)) is linear , then fora given step size it may be possible 
to find a y such that the method could be effective. If the model is nonlinear, it 
may be locally linearized to obtain a y (for a given step size) so that the method may 
remain effective. This stability analysis, done by Lomax (private communications), 
is our next topic of discussion. 


STABILITY ANALYSIS 


Let us consider a simple 


linear model : 

du 

dt 


Au 


\ 


> 


u(t ) = u I 
o o / 

If we combine equations (3) and (4) for this model we get 

U , = (1 + z + yz 2 )U 

n+1 n 

where z = AAt. Since A may be complex, so is z. It is well known that for 
stability 


( 6 ) 


(7) 


| a | < 1 

where 


(8a) 


a = 1 + z + yz 2 . (8b) 

Using z as z = x + iy , a complex variable, when we plot the equality part of 
equations (8) in a complex plane we get a stability contour for a given value of y. 
Some of these stability contours (drawn by Lomax) have been shown here in 
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the figures 1-3 for y = 0.095, 0.175, and 0.25, respectively. Outside these 
stability contours, the inequality (eqs. (8)) is not valid and as such the algorithm 
(eq. (7)) is unstable. 

If we set y = 0, we get the stability contour for the forward Euler method. 
This represents a unit circle passing through the origin and x = -2. (This is 
included in fig. 3). These contours show that knowing A, which must be real (in 
this case), and At, we can choose a y such that the method will remain stable. 
This scheme will be demonstrated in some simple examples. First we will do the 
truncation error analysis. 

In order to make a simple analysis of truncation error we consider a simple 
linear model as follows: 


. The present scheme is 


— = Au 4- a u(0) = u 

dt o 


u = u n + zu 11 4- Ata + (TE) (n) 


u n+1 = (1 - y)u + y[u n + zu + aAt 4- (TE) (n) ] (11) 

where z = AAt, (TE) = 0(At 2 ) = truncation error at t n . We know that 
lim (TE) = 0. This gives the consistency of the forward Euler scheme (eq. (10)). 
At'-O 

We will now examine how errors for the combined scheme behave as At -*• 0. 

Combining equations (10) and (11) we get 

u n+1 =. (1 + z + yz 2 )u n + (1 + yz)Ata 4- (1 4- yz)(TE)^ 

The term, neglected in the algorithm, is (1 4- y z) (TE) ^ , which tends to zero as 
At. -> 0. 

. We may study now the accuracy of the steady-state solution. Let us consider 
equation (9) again. The steady-state solution is 

• u = -a/A (13) 

Neglecting the truncation error part in equation (12) we get 

U n+1 = (1 4- z 4- yz 2 )U n + (1 + yz)Ata (14) 

where U n = the net-function corresponding to u°. 

At the steady-state condition, u” +1 = U n . Thus, from equation (14), at 


steady-state 


1T n (1 + yz)At a , . . . 

U = - , . ' \ • a = - T (z = AAt) 

z(l + yz) A 


Thus, the present numerical scheme leads to the same analytical steady-state solution 
as the one given by the differential equation itself. 
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• APPLICATIONS 

Let us consider two simple applications, one linear and the other nonlinear. 


Application No . 1 


dt^ = -100 (u 
■ u(0) = 0 


sin t) + cos t 


( 15 ) 


The analytical solution is given by 

u = sin t (16) 

Here: A = -100. If we choose At = 0.1, z = XAt = -10.0. From figure 1, z = -10 
is within the stability contour for y = 0.095. The computational results obtained 
here are given in table 1. 

Obviously, if we set y = 0,.the method reduces to the forward Euler scheme 
which was unstable for At = 0.1. There is an interesting phenomenon here: if we 

reduce the step size and choose At = 0.06, then z = XAt = -6.0 which is exterior 
to the stability contour for Y =0.095. Thus, the method must diverge. This has 
been verified computationally. 


Application No. 2: 

77 = -25 (u - 1/u) 
at ► 

u(0) = Jl 

The analytical solution is given by 

u ( t ) = [1 + exp (-50t) ] 1 / 2 

Obviously, lim u(t) = 1 which gives the steady-state solution. 


(17) 


(18) 


To obtain the value of A we may linearize the model near t = 0. Here, 
f(u,t) = -25 (u - 1/u) 

X = ~ = -25(1 +1/2) = -37.5 

du I ■ 

^-u 0 

Choosing At = 0.1, z = XAt = -3.75. Now, according to figure 1, if we choose 
y = 0.095, since z = 3.75 is outside the contour of stability, the method must 
diverge (verified computationally) . But z = -3.75 is well within the contour of 
stability for y = 0.175 (fig. 2). Computational results for this case are given in 
table 2 . 
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BURGERS' EQUATION 


One primary objective in this work is to see how effective the method could be 
with regard to the solution of mathematical models represented by nonlinear partial 
differential equations. A very common test problem for numerical methods is Burgers' 
equation, which is 


3u/3t + u 3u/3x = v 3 2 u/3x 2 . (19) 

This is a nonlinear, parabolic, partial differential equation and as such its solu- 
tion has a special significance in applied mathematics. Since, in general, closed- 
form solutions of this equation that are subject to a given set of initial-boundary 
conditions are difficult to obtain, the equation is solved numerically. 

We tested the effectiveness of our algorithm with regard to solution of 
equation (19) subject to three distinct sets of initial-boundary conditions. Some of 
our findings will be discussed now. 

Let us approximate u t by the forward Euler difference formula and the space 
derivatives u x and u xx by central differences. The explicit finite-difference 
pr edic tor is then 

U. = U n + a U n (U n , - u" ) + b(U n - 2U n + U r ‘ ) (20) 

3 3 J J-l 3 + 1 3-1 J 3 + 1 

where 

U n = U(x , t ) = the net function corresponding to u. 

3 3 n 3 

a = At/ (2Ax) , b = vAt/Ax 2 

At = time-step , Ax = net spacing 

U = predicted value of U. at t . 
j 3 n +1 

A corrector may now be set up as follows: 

u n+1 = (1 - y)U. + y[U U + a U.(U n+ ! - Cl ) + b(U n+ ’ - 2U. + li )] (21) 

3 3 3 3 3-1 .i + l 3-1 3 3 + 1 

(j = 1,2,. . . , UMAX - 1) 

where boundary conditions are given at j = 0 and j = JMAX. As before if we now 
set y = 0, then Uj + 1 = Uj , which means the corrector is not in use. We may see now 
that equation (21) retains accuracy for the steady-state solution. 

At the steady state, Uj = Uj + I = Uj for all j. Thus, if we replace in equa- 
tion (20) Uj by U 1 ?, we get the difference equation for the steady-state condition 
-• J 

as 

a U n (U n - u") + b(U n - 2U n + u" ) = 0 (22) 

3 3-1 3+1 3-i 3 3+i 

In equation (21), if we replace Uj + 1 and Uj by Uj , we get precisely the same 
equation (22). Thus, accuracy of the solution is retained by the algorithm 
(eq. (21)). 
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Let us examine now how the truncation error will behave under this convex map- 
ping analysis. Assuming that u(x,t) G C 4 ’ 2 (four times continuously differentiable 
with respect to x and twice differentiable with respect to t) we have 


u. = u° 4- G.(u n ) + (TE) n 
3 13 3 


where 


G.(u ) = a u.(u. - u., ) + b(u.. - 2u . + u. ) 

3 3 3-1 3 + 1 j + l j 3-1 

, mr ,.n . (At , .n , Ax 2 f, n v . .n"|) , 

(TE) . = At{ — (u ). + —r- (uu ) . - x- (u ) > + 

j l 2 tt 3 6 L xxx 3 2 xxxx jJJ 


(23) 


As At ■> 0 and Ax -»■ 0, (TE) -»■ 0. 


gives 


Let u .= (u ls u 2 , • • -^jmax-I^ ’ u = (u l’ u 2> 


\ U JMAX-1' > 


The corrector 


n+1 ~ , r n , „ /- n+1. , . .n, 

u = (1 - y)u. + y[u. + H.(u,u ) + (TE) . ] 

3 3 3 3 3 

Obviously (u°) may be replaced by the functional H^(u l \u 11 ) where 

H.(u,u n+1 ) = a u.(u n+1 - u ) + b(u - 2u . + u D+ j) 
3 3 3-1 3+1 3 + 1 3 3-1 


(24) 


From equations (23) and (24) 


n+ 1 n tt/ . tti/ a n+1-. . n n » , . . n 

u = u + H(u , u ) + y[H(u,u ) - H(u ,u )] + (TE) 


(25) 


If we linearize and assume that the relation 


u ,~ n+1. . n n n . n+1 n. 

H(u,u ) - H(u ,u ) = A^(u - u ) + B^(u - u ) 


(26) 


(A n and B n are square matrices of the order JMAX - 1) is valid at each time-step, 
we get, from equations (25) and (26), 


n+1 n „/ n n, , n. . n+1 n, ^ ,„.n 

u - u = H(u ,u ) + yA (u - u ) + yB (u - u ) + (TE) 

n n 

Since from equation (23) u - u n = H(u D ,u n ) + (TE) n , from equation (27) we get 

u n+1 - u n = (I - yB ) -1 (l + YA ) [H n + (TE) n ] 


(27) 


(28) 


where H n = H(u n ,u n ) . 


We assumed that (I - yBn) is invertible for all n. Since as At -+ 0 and 
Ax 0, (TE)j ->■ 0, equation (28) shows that the present scheme maintains its consis- 
tency property. It is also clear that as steady state is approached, 


(I - YB n ) -1 (I + YA n )H = 0 

giving = G^(u n ) = 0, which is true, provided (I + YA^) is nonsingular. 

In order to make some assessment of the stability properties, let us consider 
the linear Burgers' equation in the next section. An attempt will be made in 
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the future to analyze the stability properties of the nonlinear Burgers' equation, 
using explicit D-Mappings, as discussed in reference 5. 

LINEARIZED STABILITY ANALYSIS 


Let us first consider the linearized stability properties of this method. We 
consider 


3u 3u 3 2 u 

3t + C 3x " V 3x2 

Then the present method may be expressed as 


(29) 


U. = U n + a(U n - u" ) + b(U n - 2U° + U" ) (predictor) (30) 

1 J J-l 1 + 1 J-l J 1 + 1 

U n+1 = (1 - y)U . + y [U n + a(U n+1 - U. , ,) + b(U n+1 - 2U. + U. , , ) ] (corrector) (31) 
1 111-1 i + l l-i l l + l 


where 


cAt , vAt 

3 = 2^ ’ b = A^ 


(32) 


The predictor may be expressed as 


U = (b + a)U^ 1 _ 1 + (1 - 2b)ll" + (b - a)u" +1 


In matrix notation, this becomes 


U = A • U 

where A is a tridiagonal band matrix, and may be expressed as 

A = B(b + a, 1 - 2b, b - a) 

The corrector may also be expressed as 

U a+1 - y(a + b)U a ^| = (1 - y - 2by)ll. + y(b - a)lh +i + yU^ 

In matrix notation, we express this equation as 

Au n+1 = ©u + ru n 


(33) 


(34) 


(35) 


where 


= B [-y (b + a) , 

1, 0] 

(36) 

= B [0 , 1 - y - 

2by , y(b - a) ] 

(37) 

= Y I 


(38) 
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( 39 ) 


Combining equations (33) and (35) , we get 

u n+1 = a - 1 (e • a + r)u n 

The amplification matrix M is thus 

M = A - 1 (0 • A + T) (40) 

Stability is obtained if M is a convergent matrix which is true if we can find a 
norm such that 


I|m|| < i (4i) 

Some computational experiments have been conducted to find p (M) numerically; 
p (M) is the spectral radius of the amplification matrix M for some given At, Ax, 

v, and y G (0,1). (For At = Ax = 0.05, v = 10 -3 , p(M) < 1 for y = 1.01.) The 

results are given in table 3 (c = 1 in all the cases) . 

If we consider the predictor alone (eq. (30)), for stability we need b < 1/2 

and a < b. These give, respectively, 

~ 5 1/2 and ^ < 2 (c = 1) (42) 

If we consider At = Ax = 0.05, v = 0.01, the first criterion is satisfied; however, 
Ax/v = 5 > 2. 

The predictor alone is unstable, whereas the present predictor-corrector is 
stable (first row in table. 3). Also if we consider At = 0.05, Ax = 0.025, v = 0.01, 
both inequalities (eq. (42)) are violated. 

Again, whereas the predictor itself is unstable when it is combined with the 
corrector for y = 0.28 (table 3), the method becomes stable. Table 3 has been 
formed by computing p (M) , where M is given in equation (40) for a given At, Ax, 
and v, and choosing 0 < y < 1 (in general). Table 3 was done by Strate (private 
communication) . 

Fourier stability analysis will be done in the future. 


COMPUTATIONAL RESULTS 

At each t n , we first computed U-j for j = 1,2,. . ., (UMAX - 1), using the 
predictor (eq. (20)); then in a SUBROUTINE we corrected all the predicted values by 
using equation (21), where the most recent corrected values were used. 

The code and its description appeared in reference 8. It is written in 
Extended BASIC for the TRS-80 color computer made by Radio Shack by the second author. 
The boundary conditions were kept fixed at j = 0 and j = UMAX. 

Case 1 : 


u(x,0) = sin ttx 

u ( 0 , t ) = u ( 1 , t ) = 0 (43) 
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The analytical solution as 


This represents a sine wave which decays as t -»• °°. 
given in references 2 and 7 is 


where 


u(x,t) = 2ttv S 1 /S 2 

CO 

Si = n exp (-n 2 -rr 2 vt)A sin(mrx) 
n= 1 

CO 

S o = exp(-n 2 7r 2 vt)A n cos(mrx) 

n =0 

= f exp[(cos ttx - l)/2irv]dx 

- 1 

A = 2 | exp[(cos ttx - 1)/2ttv] cos(mrx)dx 
n v n 


. A 0 


(44) 

(45) 

(46) 

(47) 

(48) 


It may be noticed that for all n and v ^ 0 the integrals are convergent and 

lim Ajj = 0. The rate of convergence slows down as v gets smaller. Thus, for small 

n-*» 

values of v, it is extremely difficult to compute u(x,t) using the analytical 
expression (eq. (44)). 


The computational results obtained by the present scheme have been plotted in 
figure 4. Here v = 0.01, y = 0.25, Ax = 0.05, and At = 0.1. The time = 2.5 sec. 
Evidently, as t •> ®, u(x,t) ■* 0. (The method failed if we set y = 0.) The sharp 
overshooting near x = 1 has been caused by the boundary conditions. Such solution 
patterns are quite well known; they are not time-accurate. 

Case 2: 


u(x,0) = 1 
= 0 


0 < x < 0.1 
x > 0.1 


(49) 


This represents a moving shock. Here we chose y = 0.75 and v = 0.01. Figure 5 
describes the solution. Because of the nonconservative form of the algorithm, loca- 
tions of the shock were not traced accurately. As expected, overshooting was found 
at x = 1. It was not, however, propagated upstream. With y = 0 the oscillations 
were unbounded, which is true for Euler's forward scheme by eigenvalue analysis. In 
figure 6 an appropriate scale was used to show some of these oscillations. 


Case 3 : 


u(x,0) = <f>(x,0) 
u(0,t) = 4> ( 0 , t ) - 

u(l,t) ,= 4>(1 ,t) j 


(50) 
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where 


<j>(x,t) 


O.le A + 0 . 5e B + e C 
e _A + e~ B + e _C 


(51) 


A = (0. 5/v) (x - 0.5 + 4 . 95t) ' 

B = (2. 5/v) (x - 0.5 + 0 . 075t) • 
C = (5/ v) (x - 0.375) 


(52) 


This model is discussed in references 3 and 4. It represents the motions of two 
shocks which merge into one at the steady state. With v = 0.022 and y = 0.25, U.'s 
were computed for up to 15 time-steps and plotted in figure 7. Similar profiles ^ 
were given in references 4 and 5. With y = 0, the method failed. Oscillations are 
shown in figure 8. 

In order to remove instabilities, 3u/3x was approximated by upwind differencing 
(MacCormack) . The predictor is now 

U = U n + a U n (U n - U°) + b(U n , - 2U n + 11° ) (upwind-predictor) (53) 

3 3 33 ~ 1 3 J- 1 1 J +1 

The corrector is 

U n+1 = (1 - Y)U. + Y[U n + a U.(U n+ ) - U.) + b(U n+1 - 2U. + U. )] (54) 

3 333 J-l 3 3-1 3 3+ 1 

where a = At/Ax. Other parameters are the same as discussed before. 

Applying this formula to case 1, with Ax = 0.05, At = 0.1, v = 10 -6 , and 
y = 0.4 we obtain figure 9. Instabilities were removed; however, with y = 0, the 
method failed (fig. 10). 


The same was true for cases 2 and 3. The motion of the shock is shown in 
figure 11 for case 2, in which Ax = 0.05, At = 0.1, v = 0.01, and y = 0.25. There 
was no overshooting. However, when y = 0, the method failed after six time-steps; 
this is shown in figure 12. 

Figure 13 shows no overshooting for case 3 (Ax = 0.05, At = 0.1, v = 0.022, 
y = 0.15, number of time steps = 22). Figure 14 shows that the method collapsed 
with y = 0. 

In figure 15 it is shown that even when At = 4Ax (At = 0.2, Ax = 0.05), the 
method may still be effective for representing the motion of the shock. Figure 16 
shows somewhat interesting results for case 2. Here Ax = 0.05, At = 4Ax, v = 0.01, 
and y = 0.355. Numerical solutions initially showed some strong oscillations. 

Later they were damped out. Although the results were not time-accurate, the con- 
verged form of solution was correct. 


In figures 17-19, it is shown that for the present explicit predictor-corrector 
scheme. At, which is much larger than Ax, may be used for small values of v if 
we use the upwind predictor (eq. (27)). It has also become clear that although 
steady-state solutions are correct here, intermediary solutions at various time- 
steps are not acceptable. 
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Larger At does not necessarily mean a faster rate of convergence to steady 
state; the latter depends on how fast residuals converge to zero. This has been 
done in the past with regard to this method for solution of Burgers' equation 
expressed in the potential form. Details on this may be found in reference 9. 


DISCUSSION 


The primary objective of the predictor-corrector formulas discussed here is to 
obtain the steady-state solution. The analysis given in the Stability Analysis 
section shows that this objective is fulfilled. 

The method is only first-order accurate unless y = 1/2. Stability analysis 
for the linear Burgers' equation shows that we may use At much larger than Ax 
with a suitable y, such that the numerical solution may remain stable (table 3). 
Computationally, this has. been found to be true for the nonlinear Burgers' equation. 
In the future, nonlinear error analysis done in reference 5 using D-mappings will be 
applied to this problem. 

More analyses and computational experiments are needed in order to make a more 
realistic assessment of the present predictor-corrector schemes for numerical solu- 
tion of nonlinear models. 
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APPENDIX 


THE CODE IN TRS-80 COLOR-EXTENDED BASIC 


10 REM 

11 REM 

20 REM BURGERS' EQUATION 

25 REM 

30 REM DU/DT+U*DU/DX=NU*D2U/DX2 
35 REM 

40 REM A PREDICTOR-CORRECTOR METHOD 
45 REM 

50 REM U=VALUE OF U AT THE RECENT TIME STEP 
55 REM U0=VALUE OF U AT THE PREVIOUS TIME STEP 
57 REM G=A CONVEX PARAMETER 

60 REM JMAX=FIELD SIZE 

61 REM NTSTEPS=NO . OF TIME STEPS 

63 REM K=LEVEL OF ITERATION 

89 REM 

90 REM 

99 REM 

100 DIM U(51) ,U0(51) 

109 REM 

110 REM 

115 REM INPUT PARAMETERS 

117 REM 

120 DX=0.05 :DT=0. 1 :JMAX=21 :NTSTEPS=200 

125 NU=0. 01 :PI=3 .1416 

130 A=0.5*DT/DX:B=NU*DT/(DX*DX) 

132 K=0 

140 REM 

142 REM CHOOSE G 

145 REM 

150 G=0. 25 

160 REM 

162 REM SET INITIAL & BOUNDARY CONDITIONS 

165 REM— 

170 FOR J=1 TO JMAX :X=( J-1)*DX:U( J)=SIN(PI*X) :NEXT J 

179 REM 

180 REM 

182 REM RECORD PREVIOUS TIME STEP SOLUTIONS 

185 REM 

190 FOR .1=1 TO JMAX :U0(.T)=U( J) :NEXT J 
198 REM 

200 REM 

212 REM ITERATIONS START 

215 REM 

220 K=K+1 

229 REM 

230 REM 

232 REM EULER FORWARD : PREDICTOR 

234 REM 

240 FOR J=2 TO JMAX-1 
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245 U(J)=U0(J)+A*U0(J)*(U0(J-1)-U0(J+1))+B*(U0(J-1) 
-2*U0( J)+U0( J+l) ) 

250 NEXT J 
260 GOSUB 1000 

270 REM 

272 REM OUTPUT 

275 REM 

280 PRINT "AT TIME STEP=";K 
290 FOR J=1 TO JMAX : X= ( J-l) *DX 
295 PRINT X,U(J) 

300 NEXT J 

320 IF K<NTSTEPS THEN 190 
800 END 

900 REM 

902 REM CONVEX CORRECTOR : CHARLIE 

905 REM 

1000 FOR J=2 TO JMAX-1 
1010 TEMP=(1.-G)*U(J) 

1030 V=U0(J)+A*U(J)*(U(J-1)-U(J+1))+B*(U(J-1)-2*U(J) 
+U(J+1)) 

1050 U(J)=TEMP+G*V 
1060 NEXT J 
1090 RETURN 
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TABLE 1.- COMPARISON BETWEEN 
ANALYTICAL AND COMPUTATIONAL 
VALUES OF U: v = 0.095, At = 0.1 


Time 

U (Analytical) 

U (Computed) 

2.5 

0.598472144 

0.597493063 

5.0 

-0.958924274 

-0.957553061 

7.5 

0.937999977 

0.936781981 

10.0 

-0.544021111 

-0.543440746 

12.5 

-0.0663218955 

-0.0660338101 

15.0 

0.650287837 

0.649245878 


TABLE 2.- COMPARISON BETWEEN 
ANALYTICAL AND COMPUTATIONAL 
VALUES OF U: v = 0.175, At = 0.1 


Time 

U (Analytical) 

U (Computed) 

0.1 

1.00336332 

1.12706991 

0.5 

1.00 

1.04808666 

0.9 

1.00 

1.00621156 

1.3 

1.00 

1.0001474 

1.7 

1.00 

1.00000293 

2.1 

1.00 

1.00000006 
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TABLE 3.- STABILITY FOR LINEAR BURGERS' 
EQUATION 



c = 

1.0 

v = 0.01 

Case 

At 

Ax 

P MIN 

GAMMA 

n 

1 


0.6741081217 

ESI 

IS 


0.025 

0.2634833436 

■ 

11 


0.05 

0.5591445292 

m ■ 

n 


0.025 

0.4932129687 

0.22 

■ 

0.2 

0.05 

0.8777235062 

0.32 

H 

0.2 

0.025 

0.731106526 

0.13 


c = 

1.0 

v = 0. 

L 

i 

0.05 

0.05 

0.8349436353 

0.18 

2 

0.05 

0.025 

3.8000485810 

0.05 

3 

0.1 

0.05 

1.5907523689 

0.1 

4 

0.1 

0.025 

8.0442946351 

0.03 

5 

0.2 

0.05 

3.7391180937 

0.05 

6 

0.2 

0.025 

24.990423418 

0.01 


c = 

1.0 

v = 0.001 

1 

0.05 

0.05 

0.95710316836 

1.010 

2 

0.05 

0.025 

0.8299183481 

0.77 

3 


0.05 

0.8320843664 

0.89 

4 


0.025 

1.7710418247 

0.37 

5 

0.2 

0.05 

1.9420702213 

0.38 

6 


0.025 

4.5407595319 

0.18 
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CONVEX CORRECTOR 7 = 0.095 



Figure 1.- Stability contours, convex corrector: y = 0.095. 


CONVEX CORRECTOR 7 = 0.175 
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Figure 19.- Velocity profiles using upwind differencing for case 1: v = 0.00001, 

y = 0.152, Ax = 0.05, At = 0.3, number of time steps = 40. 
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